function g = maxgrad(d)
% maxgrad  - return maximum direct neighbor gradient (diff) per voxel
%
% FORMAT:       mg = maxgrad(d)
%
% Input fields:
%
%       d           3D data
%
% Output fields:
%
%       mg          max. gradient

% Version:  v0.7f
% Build:    8110521
% Date:     Nov-05 2008, 9:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin ~= 1 || ...
    ndims(d) ~= 3 || ...
   ~isnumeric(d)
    error( ...
        'BVQXtools:BadArgument', ...
        'Bad or missing argument in call.' ...
    );
end

% make double
if ~isa(d, 'double')
    d = double(d);
end

% create three convolution operators
cop = zeros(3, 3, 3);
cop(2, 2, 2) =  1;

% compute first direct gradient
cop(1, 2, 2) = -1;
g = abs(conv3d(d, cop, 2));
cop(1, 2, 2) = 0;

% then max with others
cop(3, 2, 2) = -1;
g = max(g, abs(conv3d(d, cop, 2)));
cop(3, 2, 2) = 0;
cop(2, 1, 2) = -1;
g = max(g, abs(conv3d(d, cop, 2)));
cop(2, 1, 2) = 0;
cop(2, 3, 2) = -1;
g = max(g, abs(conv3d(d, cop, 2)));
cop(2, 3, 2) = 0;
cop(2, 2, 1) = -1;
g = max(g, abs(conv3d(d, cop, 2)));
cop(2, 2, 1) = 0;
cop(2, 2, 3) = -1;
g = max(g, abs(conv3d(d, cop, 2)));
